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Large-scale non-Gaussian mass function and halo bias: 
tests on N-body simulations 
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ABSTRACT 

The description of the abundance and clustering of halos for non-Gaussian initial con- 
ditions has recently received renewed interest, motivated by the forthcoming large 
galaxy and cluster surveys, which can potentially yield constraints of order unity on 
the non-Gaussianity parameter Jnl- We present tests on N-body simulations of ana- 
lytical formulae describing the halo abundance and clustering for non -Gaussian initial 
condit ions. We calibrate the anal ytic non-Gaussian mass function of iMatarrese et all 
()2000( ) and iLoVerde et al.l (|2008l ) and the analytic description of clustering of halos 
for non-Gaussian initial conditions on N-body simulations. We find excellent agree- 
ment between the simulations and the analytic predictions if we make the corrections 
Sc — > Sc X y/g and Sc — > Sc x q where q ~ 0.75, in the density threshold for 
gravitational collapse and in the non-Gaussian fractional correction to the halo bias, 
respectively. We discuss the implications of this correction on present and forecasted 
primordial non-Gaussianity constraints. We confirm that the non-Gaussian halo bias 
offers a robust and highly competitive test of primordial non-Gaussianity. 

Key words: methods: statistical, iV-body simulations-cosmology: theory, large-scale 
structure of universe galaxies: clusters: general - galaxies: halos 



1 INTRODUCTION 

Constraining primordial non-Gaussianity offers a power- 
ful test of the generation mechanism of cosmological per- 
turbations in the early universe. While standard single- 
field models of slow-roll inflation lead to small departures 
from Gaussianity, non-sta ndard scenario s allow for a larger 
level of non-Gaussianity iBartolo et al.l ([2004) and refer- 
ences therein). The standard observables to constrain non- 
Gaussianity are the Cosmic Microwave Background (CMB) 
and the Large Scale Structure (LSS) of the Univers e. A pow - 



servational aspects of non-Gaussianity have been exten- 
sively tested against N-body simulation s, leading to dif 



ferent and sometimes conflicting results (iKang et al 



Grossi et aTlbOOTl: lOalal et al.|[2007l : lDesiacaues et al 
Pillepich et al.ll2008l V 



2007 



2008 



erful technique is based on the abunda nce ( Matarr ese et al 
2000| : IVerde et al.ll 2001: Lo Verde et a l. 2008; Kova ma et al 
1999 : iRobinson fc Baker .2005 



Robmson et al. 20001) and 



- ___oir 

clustering (ICrinstein fc Wise! Il986l : IMatarrese et al.l [19861 : 



iLucchin et al.l 119881 ) o~ rare events, such as dark matter 



density peaks, as they trace the tail of the underlying 
matter distribution. Theoretical predictions on various ob- 



iDalal et al.l (|2007l ) and IMatarrese fc Verdd (|2008l ) 

showed that primordial non-Gaussianity affects the cluster- 
ing of dark matter halos inducing a scale-dependent bias 
on large scales. Not only this effect has been already ex- 
ploit ed to place st r ingent constraints on no n-Gaussianity 
(Slos ar et al] l2008l : lAfshordi fc ToUevI l2008f ). but also it 
is particularly promising for constraining non-Gaussianity 
from future surveys, which will provide a large sample of 
galaxy clusters over a volume comparable to the horizon size 
(e.g., DES, PanSTARRS, BOSS, LSST, ADEPT, EUCLID) 
(Dalai ct al. 2007; Carbonc ct al. 200a lAfshordi fc ToUevI 
I2OO8I : ISehakH2008i '). iBartolo et all l|2005l ) showed that even 



© 0000 RAS 



2 Grossi et al. 



for small primordial non-Gaussianities, the evolution of per- 
turbations on super-Hubble scales yields extra contribu- 
tions. The amplitude of these contributions is comparable 
to the forecasted errors of some planned surveys, opening 
up the possibility of measuring them. 

In light of this, it is important to use N-body exper- 
iments to test the validity of theoretical predictions for 
halo-bias in non-Gaussian framework. Indeed, all proposed 
analytic biasing expressions have been derived in the ex- 
tended Press-Schechter framework which assumes spherical 
collapse dynamics, sharp k-space filtering and Gaussian ini- 
tial conditions. The validity of the extrapolation of the ex- 
tended Press-Schechter approach to the non-Gaussian case 
can be tested independently by considering also the halo 
mass function. It is thus also important to test and cal- 
ibrate on N-body simulation s the pr edictions of the non- 
Gaussian halo mass fun ction (|Kang et a l. 2007.: lGrossi et al.l 
I2OO7I : iDalal et al.ll2007l ) and of the non-Gaussian halo bias 
simultaneously. This is what we set out to do here. 

In this paper we start by reviewing the analytic pre- 
dictions for the Gaussian and non-Gaussian halo abundance 
and clustering (§2). In §3 we describe the numerical simula- 
tions with Gaussian and non-Gaussian initial conditions. In 
§4 we present the test for the non-Gaussian mass function. 
In §5 and 6 we test the analytic predictions of Gaussian and 
non- Gaussian large scale bias against N-body simulations. 
In §7 we compare our results with the literature. Finally, we 
conclude in 58. 



2 FORMULATION OF THE NON-GAUSSIAN 
HALO ABUNDANCE AND CLUSTERING 

Deviations from Gaussian initial conditions are commonly 
par ameterized in ter ms of the dimensionless /nl parame- 
ter JSalopek fc Bo nd 1990; Gan gui et aT]|l994l : IVerde et all 



I2OOOI : iKomatsu fc Spergel.2001 ): 



$ = <^ + /nl(0'-{0')), 



(1) 



where "I> denotes the gravitational po tential and </) is a Gaus - 



sian random field. As no ted by e.g., iLoVerde et al.l 1 200; 



lAfshordi fc ToUevI l|2008D and iPillepich et al.1 (|2008l ). dif- 



ferent authors use different conventions. Here $ denotes 
Bardeen's gauge-invariant potential which, on sub-Hubble 
scales, reduces to the usual Newtonian peculiar gravitational 
potential but with a negative sign. In addition, there are 
two conventions for normalizing Eq. the LSS and the 
CMB one. In the LSS convention $ is linearly extrapo- 
lated at 2 = 0. In the present paper we use this conven- 
tion. In the CMB convention "1> is instead primordial: thus 
/nl = [g{z = oo)/s(0)]/^r'' ~ ^, where g{z) de- 

notes the linear growth suppression factor in non Einstein- 
de-Sitter Universes. 



2.1 Formulation of the non-Gaussian mass 

function: Extended Press-Schechter approach 

In the Press-Schechter framework, one considers the density 
contrast field evaluated at some early time, far before any 
scale of interest has approached the nonlinear regime, but 
extrapolated to the present day using linear perturbation 
theory. Then one considers the height of the critical density 



threshold as a function of time. In that way, the collapse of 
a halo at redshift z ^ corresponds to the z = density 
fiuctuation crossing a barrier of height Sc{z) — AcD{z = 
0)/D{z), where Ac ~ Sc{z — 0) (this is an equality only 
in an Einstein de Sitter Universe); we use D{z = 0) = 1, 
D{z) = g{z) / g{0){l + z)~^ . We should recall here that, even 
in linear theory, the normalized skewness of the density field, 
S3 = {S^)/{5^)^, depends on redshift cx 1/D{z), however 
in the Press-Schechter framework one should use the linear 
5*3(2 = 0), in what follows S3 = S-s{z — 0). Note also that 
in general the skewness can be written as S3 = /nlSj^', 
where Sg^' denotes the skewness in units of /nl, care must 
be exercised in the interpretation of /nl: if Sg^' is that of 
the density field linearly extrapolated at 2; = 0, /jvl must 
be the LSS one and not the CMB one. 

Generalization of the mass function to non-Gaussian 
initial conditions within the Press-Shechte r formalism has 
been presented in iMatarrese et al.| l|2000l '): iLoVerde et al.l 
(|2008l ). Both references start by computing an expression 
for the non-Gaussian probability density function of the 
smoothed dark matter density field, then obtain the level 
excursion probability. In the Press-Shechter approach the 
mass derivative of the level excursion probability is the key 
ingredient to obtain the mass function expression and is the 
term that gets modified in the presence of primordial non- 
Gaussianity. In this derivation, several approximations are 
made. Both approaches assume that deviations from Gaus- 
sian ity are small. 

IMatarrese et al.l l|2000l ) use first the saddle-point ap- 
proximation to compute the level excursion probability and 
then truncate the resulting expression at the skewness. They 
obtain Q: 



n{M,z) = 2 



m,0 



exp 



dS3,. 



6', 



(2) 



+ 



5f do" A 



2 3^1-S3,A/5c/3'^lnM audlnM 

where gm denotes the rms value of the density field, the sub- 
script M denotes that the density field has been smoothed 
on a scale R{M) corresponding to R{M) = [M3/(4pA/)]^''^, 
and St = Scy/l — 5 c S3,m/ 3. 

ILoVerde et al.l l|2008l ) instead first approximate the 
probability density function using the Edgeworth expansion, 
then perform the integral of the level excursion probability 
exactly on the first few terms of the expansion. They obtain: 



n{M,z) = 2 



exp 



2al 



(3) 



d In aM / _5c_ , Sz^MO-M 



dM 



Vcrjvi 



- + 



, 1 dSs^M 



Note that in the limit of small non-Gaussianity and rare 
events, the ratio of the non-Gaussian mass function to the 
Gaussian one for both expressions reduces to: 



1 We correct here a typo in Eq. (68) of IMatarrese et al] ll200(t) . 
where diner ]\i should be dcrj^f. 
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TZng 



n(M, 2|/n 



1 + Ss 



(4) 



n(M, 2|/nl = 0) 

It is important to bear in mind that in Eqs. ©-Q, the 
redshift dependence is enclosed only in 5c (and not in 5*3). 
In the spirit of the "CMB" convention instead, where the 
gravitational potential is normalized deep in the matter era, 
one should make sure that all the relevant quantities are cor- 
rectly extrapolated linearly at z = 0, keeping in mind that 
the gravitational potential slowly evolves in a non Einstein 
de Sitter Universe. 

The major limitations in both derivations are the as- 
sumption of spherical collapse and the sharp fc-space filter- 
ing. In addition, the excursion set improvement on the in- 
terpretation of the original Press-Shechter swindle, suggests 
that this de rivation relies on the random-phase hypothesis 
IShethI I 19981 ) ■ which is clearly not satisfied for non-Gaussian 
initi al conditions even for s harp fc-space filtering . 

IVerde et al] ([2001) and lLoVerde et al.1 (|2008l ) addressed 
this issue by using the analytical approach to compute the 
fractional non- Gaussian cor rection to the Gaussian mass 
function TZng, and used the ISheth fc TormenI l|l999D mass 
function to model the Gaussian mass function. This ap- 
proach is potentially promising, but needs to be calibrated 
on numerical experiments. 

In particular, one may argue that the same correction 
that in the Gaussian case modifies the coll apse threshold and 
thus t he f orm of the m ass function from Press & Schechter 
l| 19741 ) to ISheth et al] (2001) and Sheth & Tormcn (2002), 
may apply to the non-Gaussian correction. In the Gaussian 
case this is usu ally referred to as the c orrection due to ellip- 
soidal collapse (|Lee fc Shandarinll 19981 '). W hile this interpre- 
tation has recently been disputed (see e.g.. [Robertson et al.l 
l|2008l )l, we will maintain the same nomenclature here. For 
rare events, high peaks {5c/ aM 3> 1) and small /nl, this is 
equivalent to lower 5c by a factor with q = 0.75. 

In summary we propose that the non-Gaussian mass 
function n{M, z, /nl) should be re- written in terms of the 
Gaussian on e nfi'^(M,z) -gi v en by t ested fits t o simu - 
lations e.g.. 



[Sheth fc TormenI lll999[): iReed et all l|2003D : 

IWarren et all (120061 ): I Jenkins et al.l l|200lh -, muhiplied by 
a non-Gaussian correction factor: 



n{M,z,fML) = n(i'"(M,z) x TlNoiM, z, /nl) 



(5) 



where Rng(M, z , /nl ) tak es two differen t form s in the 
iMatarrese et al.1 hOO&i and 'LoVerde et alj (I2OO8I ) approx- 
imations. For the [Matarroso ct al. (2000) case[f| 



RNGiM,Z,fNL) 

sic 



1 



exp 



dS-i 



Sz.M 



(6) 



5ec\/l — 5ecS3,M/3 



5ccSi,M/'id\r^aM 
and for the lLoVerde et aL I l|2008t ) case: 



Rng{M,zJnl) 



st, 



1 + i^x 

6 5ec 



(7) 



Sz.M [^-2^-1 



dSz.M 
d In (jm 



- 1 



We correct here a typo in Eq. (3) of lGrossi et al. I iliooj) where 
the exponential part was missing. 



where 5cc denotes the critical density for ellipsoidal collapse, 
which for high peaks is 5cc ~ 5c-^ with q — 0.75. 



2.2 Formulation of the non-Gaussian large scale 
halo bias 

For the case of "local" primordial non-Gaussianity Eq. (|lj, 
the analytical expression for the large-scale non-Gaussian 
bias has been derived in fi ve different ways, o btaining always 
basically the same resuh. iDalal et al.1 (|2007l ) considered the 
Laplacian of $ in the vicinity of rare, high peaks, consider- 
ing that the resulting is proportional to the peaks over- 
density; they also generalized to local non-Gaussianity the 
iKaiseil l|l984 ) argument o f high-peaks bias in order to derive 
its non-Gaussian version. IMatarrese fc Verd3 l|2008l ) derived 
the halo bias formula in general non-G aussian cases speci- 
fied by an expression for the bispectrum. lSlosar et al.l (l2008f) 
adop ted the peak-background split approach ( Cole fc Kaiserl 
Il989l ) for the local non- Gaussian case, showing that the 
resulting expression relies on the universality of the mass 
function. lAfshordi fc TollevI ( |2008| ) instead interpreted non- 
Gaussianity as a modification of the critical density for 
collapse, in the fr amework of ellipsoidal collapse. Finally, 
iMcDonald I l|2008l ) used a renormalized perturbation the- 
ory approach to consider at the same time non-linear 
bias, second-order gravitational evolution and local form of 
non-Gaussianity. It is encouraging that these different ap- 
proaches yield a consistent result for the correction to the 
Gaussian Lagrangian halo bias 6^: 

A6 

W 



2/NL(5c(2)Q!Af (fc) 



(8) 



where amik) encloses the scal e and halo mass dependen ce 
-see e.g., Eq. (13) and Fig. 3 of lMatarrese fc Verdel (|2008l )-. 
Also in this case the density field is the one extrapolated 
linearly at z = 0, and am does not depend on redshift. 

Making the standard assumption that halos move co- 
herently with the underlying dark matter, the Lagrangian 
bias is related to the Eulerian one as & = 1 4- . 

The approxi mations used to derive this equation are 
Press-Schechter (|Press fc Schechter! Il974l ) approach, linear 
bias, small non-Gaussianity, and in most cases spherical col- 
lapse and identification of peaks with halos. It is therefore 
important to test the validity of Eq. (jSjl, with simulations 
and see if any correction factor needed is indeed due to ac- 
co unt for non-spheric a l coll apse. Following the derivation 
of lMatarrese fc Verdel l|2008l ) we recognize that the correc- 
tion to the 2-point halo correlation function due to non- 
Gaussianity (their Eq. (6)) is multiplied by v'"^ /(j'm with 
!^ = 5c/o"Af . In this factor we recognize one Lagrangian Gaus- 
sian bias factor to the second power and an extra 5c/o-\i, 
which denominator was absorbed in the form factor. Re- 
call that, as discussed in §2.1, for "ellipsoidal collapse" and 
rare events, the Lagrangian Gaussian bias is corrected as 
V jcFM — > qv/cTM (see Eq. pip below, for high v). However, 
the remaining factor is also a Gaussian bias and it should 
also be corrected by the g-factor. 

We conclude that the "non-spherical collapse" modifies 
Eq. dH} to be: 



^ ~ 2/NL<5c(2)QAf(fc)? . 
^1 



(9) 
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Figure 1. Skewness S3 of the smoothed initial density field for 
/nl = 100, 200, 500, 1000. Symbols show the numerical results 
of the initial conditions code (averaged over 5 realizations) and 
are plotted against the analytical predictions for smoothing radii 
Ts = 4, 6, 8, lOMpc/h of a spherical top-hat filter. 



10 



0.01 < k< 0.03 



V=0o 

fNL=^OOA 

C=-ioo 

C=-200 



10 



Figure 2. The bias of the halo power spectrum 6^^ compared 
to the bias of the cross (halo-matter) power spectrum bhm- As 
expected, when the number density of halos is high there is good 
agreement between the two quantities. At low halo number densi- 
ties the two quantities are affected differently by shot noise, with 
bhm being the least affected. 



Note that lAfshordi fc TollevI | 20081 ) arrived to a similar yet 
not identical expression when considering ellipsoidal col- 
lapse, i.e. they sugg est that dc s hould be substituted by 
the critical density of lSheth et al.l l|200ll ). which in our limit 
would correspond to use rather than q in Eq. Q . 

In SjS]we will show that Eq. (|5J correction fits well the 
simulations. 



3 N-BODY SIMULATIONS 

The deviations from Gaussianity we are after become im- 
portant on very large scales k < 0.03h/Mpc and for massive 
halos. Therefore, one needs to perform N-body simulations 



on very large boxes, yet with enough resolution to identify 
massive virialized structures at different redshifts. 

Suitable initial conditions have be en set up fo ll owing 
the meth od described i n mor e detail in Grossi et al.l (|2008l ) 
(see also iGrossi et al] (|2007l ): IViel et al.l (|2009l )). In brief, 
a random realization of a Gaussian gravitational potential, 
<E>L, normalized to be the one linearly extrapolated ai z — 0, 
is generated in Fourier space, then it is inverse-Fourier trans- 
formed back to real space and added to the non-Gaussian 
term, <1>nl = /nl ("I-l - (^l)) • The resulting field $l -|-$nl 
that is linear and at 2; = 0, is transformed back in Fourier 
space. We eventually modulate the power-law spectrum us- 
ing the transfer function and compute the corresponding 
density field, which we then scale back to the initial con- 
ditions redshift {z — 60). The corresponding gravitational 
potential is then used to displace particles according to the 
Zel'dovich approximation. This method allows one to sim- 
ulate non-Gaussian models having power spectra which are 
all consi stent with t hat of the Gaussian case and was already 
used by lViel et all (|2009l ). 

In order to check the reliability of the initial condi- 
tions generation, we have performed a specific test: using 
256"^ particles in a box of size 1000 Mpc/h, primordial den- 
sity fields (extrapolated linearly at z = 0) were generated 
and smoothed using spherical top-hat filters of different radii 
Vs = 4, 6, 8, 12 Mpc/h. The smoothed skewness was then 
extracted from the fields and compared to the analytical 
prediction for /nl ~ 100, 200, 500, 1000, as shown in Figure 

m 

The set of simulations used in this work assumes the 
'concordance' ACDM model. We fix the relevant parame- 
ters consistently with those derived from t he analysis of the 
WMAP 5-year data (|Komatsu et al.ll2008t ): Qmo = 0.26 for 
the matter density parameter, Hao = 0.74 for the A con- 
tribution to the density parameter, h — 0.72 for the Hub- 
ble parameter (in units of 100 km s^^ Mpc~^). The ini- 
tial power spectrum adopts t he Cold Dark Matter (C DM) 
transfer function suggested by lEisenstein fc Hul l| 19991 ) . has 
a spectral index n = 0.96 and is normalized in such a 
way that erg = 0.8. In all expe riments, perfo rmed using 
the GADGET-2 numerical code (ISpringelll20"05l ). switching 
off the hydrodynamical part, we consider a box of (1200 
Mpc/Zi)'' with 960"^ particles: the corresponding particle 
mass is then m ~ 1.4 x 10^^ Mq . The gravitational force 
has a Plummer-equivalent softening length of e; = 25h~^ 
kpc. The runs produced 15 outputs from the initial redshift 
{z = 60) to the present time. The 5 simulations consider dif- 
ferent amounts of primordial non- Gaussianity, parametrized 
by the /nl parameter: /nl = (i.e. the reference Gaussian 
case) and /nl = ±100, ±200. The catalogues of dark matter 
haloes are extracted from the simulations using the standard 
friends-of- friends algorithm adopting a linking length of 0.2 
times the mean interparticle distance; only objects with at 
least 32 particles are considered. 

We thus measure the halo bias in the simulations as 



bs{k,M,z) = 6h„ 



Phm{k,Z,M) 



(10) 



where Phm{k,z,M) denotes the cross-power spectrum of 
dark matter with halos of mass M at scale k, and for the 
simulation snapshot at redshift z. Similarly Pmm{k,z) de- 
notes the dark matter power spectrum. Here and hereafter 
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Figure 3. Multiplicity mass function for the Gaussian simu- 
lation computed using a Friends-of-Friends halo finder. Points 
denote the simulations results at different redshift: z = 
0,0.44,1.02,1.53,2.26 and 3.23 (top to bottom). Solid (green) 
lines are the Sheth & Tormen (1999) formula, dashed (red) lines 
are the Warren et al. (2006) one and dotted (blue) are the Jenkins 
et al. (2001). 



the subscript s denotes quantities measured from the simu- 
lation. 

In principle, the quantity one is interest ed in woul d be 
the bias of the halo power spectrum b^h ~ \/ Phh /Pmm , but 
hhm is a less noisy quantity (the shot noise of the finite num- 
ber of halos is greatly suppressed in the estimate of Phm{k)). 
The quantity bhm is not guaranteed to be identical to bhh 
if bias has a stochastic component that does not correlate 
with the matter density field. In Fig. [2] we show that this 
is not the case and that there is good agreement on large 
scales between b^h and bhm, justifying using the less noisy 
bhm as an estimator for bhh- 



3.1 Comparison with independent simulations 

In Fig. |3] we show the mass function extracted from 
our Gaussian simulations at the following redshifts: z = 
0.0,0.44,1.02,1.53,2.26 and 3.23. We also show three dif- 
ferent theor etical predictions (also c alibrated on N-bod y 



simulations) 



IWarren et al 



Sheth fc Tormenl(|l999l ): I Jenkins et al.1 (|200ll ): 
120061). solid, dotted and dashed lines respec- 



tively. There is good agreement even at high redshift. 

Several groups recently presented N-body simulations, 
aiming at quantifying the effect of the non Gaussian ini- 
tial conditions on the halo mass function l|Desiacaues et ahl 
l2008l : iPillepich et~al] |2008| ). All these results are obtained 
for similar cosmological parameters, so that we can compare 
estimates derived from all the simulations directly. By com- 
paring the results for the individual simulations at 2 = 1, 
z ^ 0.5 and z = in Figures |4l [5] we demonstrate that these 
results are in agreement among the difi'erent groups, once the 
/jvL values are suitably converted to the same convention. 
Although all simulations use boxes of Giga parsec scales to 
explore the effect of non-Gaussian initial conditions at the 
high mass end, the statistical errors at the scale of massive 
clusters are still large. Therefore, we also report the recip- 



Figure 4. Comparison between the halo mass function recov ered 
in our simulati o ns wi th the work of iDesiacgues et al. I 1 I2OO8I) and 
IPillepich et al. at 2: = 0. We show the ratio between our 

non-Gaussian and Gaussian simu lation with fisiL = ±10 0, few 
points we read out from Figure 1 oflPesiac gues et al.lll2008 l ) (black 
points) at the values of u corresponding to 1 X 10^^,1 X 10^ ^ and 
1 X IOI5 Mq/H and the points from IPillepich et al.l 1I2OO8I ). We 
plot the reciprocal of the results for Jml = —100. 



rocal of the results obtained for negative /nl so that they 
appear in the positive part of the plot, to give an intuitive 
feeling of the noise within the individual simulations. 

In Figure |4j we show our simulation results for /jvl — 
100 (blue triangles) and for /jvl = —100 (red squares) 
at 2 = comp a red w ith data points from Figure 1 of 
jPesiacques et al] (I2OO8I ) (black points) at the values of u 
corresponding to 1 x 10^^,1 x lO" and 1 x lO^'^ Mp/fc (as 
given in their figure caption). Note that, as lDesiacgues et al.l 
(|2008l ) use /jvL = 100 in the CMB convention for their sim- 
ulations, we scaled the points down accordingly by a factor 
1.3 to be com parable with ou r /nl = 100. We also show 
the results for iPillepich et al.l l|2008h (green points) . Here 
we again apply the re-scaling as before, as their /jvl of 82 
would correspond to a /at z, of ~ 106 in the LSS notation. 

In Figure [5l th e left panel shows the results for 
jPillepich et al.l (|2008l ) (green points) at z = 0.5, and our 
points for the two closest available output times of our sim- 
ulation (2 = 0.44 and z — 0.61). The right panel shows the 
comparison at 2 = 1 betwee n our points (blue triangles and 
red squares) and points from lDesiacgues et al.l |2008f ) (black 
squares) . 

/.From this comparison we conclude that there is re- 
markable agreement between the three independents simula- 
tions, highlighting the robustness of the simulations results. 
The differences visible at some of the highest mass bins are 
not significant, given the large error bars present. 



4 MASS FUNCTION 

We compare the halo mass function of the non-Gaussian 
simulations with the theoretical predictions of Eqs. ((5)1, ([6| 
and (O that is, including our ansatz for the the non-spherical 
collapse correction: 5c — > y^^c- For clarity we show here 
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Figure 5. Co mparison between the halo mass function recovered in our simulations with the work of iDesiacgues et al. I 1 I2OO8I) and 
iPillepich et al.l 12008). In the left panel we show the ratio between the non-Gaussian and Gaussian simulation at redshift 1 for our 
simulations with /jvl = ±100, three points we read out from Figure 1 of IDesiacgues et al.l ll2008h (black points) at the values of u 
corresponding to 1 X 10^'',1 X 10^* and 1 X 10^^ Mq/H. We plot the reciprocal of the results for /nl = —100. In the right panel we show 
the data of Pillepich, Porciani, Hahn (2008 ) at 2 = 0.5 and we com pare them with our simulation results for the two closest available 
redshifts : 2 = 0.44 and z = 0.61 and with IDesiacgues et al ] l|2008l V All points are rescaled to I/atl = 100| in our notation. The three 
independent simulations are in good agreement. 



the non-Gaussian to Gaussian mass function ratio, i.e. the 
factor Rmg{M, z). The comparison between theory and sim- 
ulations results is shown in Fig. |S] for a few redshift snap- 
shots and for /jvi = ±100, and in Fig. [7] for /jvi = ±200 
for the same red s hifts. Dashed lines are the mass function of 
iMatarrese et al.l (l200(]|)-Eq. ^ - and dot-dashed lines are 



2p 



that of ILo Verde et al.r(l2008l) - -Eq. Q - 

Contrary to lKang et all (|2007|) andjPalal et"ai] l|2007l ). 



we conclude that both IMatarrese et al. I (|2000l ) and 
ILo Verde et all (|2008l ) are good descriptions of the non- 
Gaussian correction to the mass function, once the correc- 
tion for non-spherical collapse is inclu ded. 

Figs.[6]and[7]seem to indicate that lLoVerde et all (120081) 
may b e a better fit for small masses and Matarrese et al.l 
l|200(]| ) at high masses. This is not surprising: the Edge- 
worth expansion works well away from the extreme tails of 
the distribution (i.e. for mod erate Sc/ctm), wh i le the saddle- 
point-approximation used in IMatarrese et al.| l|2000l ) , is ex- 
pected to work better at the very tails of the distribution 
(very high Sc/cfm)- W e expect that the mass function of 
IMatarrese et al.l l|2000l ) will be a better fit at very high 
masses or larger /jvl. This will be further explored in fu- 
ture work. 



5 GAUSSIAN HALO BIAS, AND THE EFFECT 
OF MERGERS 

The lar ge-scale, linear halo Eulerian bias for th e Gaussian 
case is ||Mo et al.ll 19971 : IScoccimarro et al.|[200ll ) 



1 + 



D{zo) 



qScizf] 



5c{zf) 



+ 



(11) 



ScMDiZo) 



qSlizf) 



where q = 0.75 and p — 0.3, account for non-spherical 
coUapse and are a fit to numerical simulations. Here, aM 
denotes the rms value of the dark matter fluctuation fleld 
smoothed on a scale R corresponding to the Lagrangian ra- 
dius of the halos of mass M; Zf denotes the halo formation 
redshift and Zo denotes the halo observation redshift. As we 
are interested in massive halos, we expect that Zf ~ Zo- 
As the non-Gaussian halo bias correction is proportional to 
— 1, the dependence of b'^ on whether the selected ha- 
los underwent a recent merger (i.e. Zf ~ Zo) or are old halos 
(i.e. Zf Zq) affects the amplitude of the non- G aussian cor- 
rection (jSlosar et al.ll2008l : [Carbone et al.ll2008l ). Before we 
trust our simulation to accurately describe the non-Gaussian 
halo bias we check whether we recover the Gaussian one 
and whether the linear halo bias approximation is a good 
description f or the scale s , reds hifts and mass-ranges we are 
interested in. iGao et all (2OO5I) show that analytical predic- 
tions for the Gaussian halo bias are in reasonable agreement 
with simulations and that the bias for low-mass halos shows 
strong dependence on formation time but high mass halos 
(the ones we are interested in) do not. The halo bias for the 
Gaussian simulation and the comparison with the theory 
prediction is shown in Fig. |8] Except for the Gaussian halo 
bias bp = l + 5c(zo) / (a%rp (zo)) defined in lEfstathiou et alTI 
(Il988t ) and iKaiseil (|l984l ) indicated by the dotted (blue) 
line, the simulated data agree with the theoretical expecta- 
tions at different redshifts. In particular, in Fig.|8l the black 
solid line represents the total Gaussian bias of Eq. (|lip . the 
dashed (red) line represents the contribution from the first 
line of Eq. (|ll|l . and, finally, the dot-dashed (green) line is 
1 + qi^o ~ 1)- The small difference when using Zf ~ Zo im- 
plies that, for the Gaussian halo bias of very massive halos 
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Figure 6. Ratio of the non-Gaussian (/jvz, = ±100) to Gaussian mass function for different redshift snapshots: top left z = 0.61; top 
right z = 1.02; bottom left z = 1.53; bottom right z = 1.86. The dashed line is the mass function of Matarrese, Verde &; Jimenez (2001) 
and the dot-dashed lines are that of Lo Verde et al. (2008), both including the g-correction. 



(M > 10" Mq), it is reliable to assume that the correction 
from the "non-spherical collapse" can be encapsulated in the 
factor q in front of Sc{zo) / {um D{zo)) ■ 



6 NON-GAUSSIAN HALO BIAS 

Eq. (ISJ shows that the redshift and scale dependence of the 
non- Gaussian correction can be factorized as a term that 
depends only on redshift and one that depends only on k 
and M. The M-dependence is expected to be very weak at 
large scales {k < 0.03 h/Mpc). Here we will test the mass, 
scale and redshift dependence of the non-Gaussian halo bias 
and we calibrate its normalization on the simulations. 

In Fig. |9] we show the dependence on halo mass of 
Ab/bL- We define the quantity 

\ Ol / s \ Ol / theory 

where {Ab/bL) theory is given by Eq. (|8}. To study the 
mass dependence, we evaluate the theory at fixed mass 
M — 10^'^ Mq. We compute the bias from the simulations 



14 - 

12- bo 

- 

10 r (b°-1)q+1 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 



Figure 8. Black solid line: the large-scale Gaussian halo Eulerian 
bias b'^ of Eq. {TlJ. Blue dotted line: the Gaussian halo bias = 
1 + Se{zo)/{'^]^j D{zo)) as defined in Efsthatiou et al. (1988) and 
Kaiser (1984). Green dotted-dashed line: l+q{b^ — 1). Red dashed 
line: the contribution bp = 1 + [qSc{z f) / aj^j — 1/Se{zf)]/ D(zo) to 
the total bias of Eq. Ulljl . 
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Figure 7. Ratio of the non-Gaussian (/jvz, = ±200) to Gaussian mass function for different redshift snapshots: top left z = 0.61; top 
right z = 1.02; bottom left z = 1.53; bottom right z = 1.86. The dashed line is the mass function of Matarrese, Verde &; Jimenez (2001) 
and the dot-dashed lines are that of Lo Verde et al. (2008), both including our q-correction. 



taking halos in six diflterent mass bins. Fig. [5] includes only 
scales k < 0.03 ft/Mpc, different lines correspond to different 
redshift snapshots between 2 = and z — 1.5. As expected, 
there is no noticeable dependence on halo mass. 

Having confirmed the expected weak dependence on 
halo mass for masses M > IO^^Mq/H and on scales k < 
0.03 ft/Mpc, we can study the redshift and scale depen- 
dence of Ab/bL, considering halos of different masses above 
IO^^Mq/Zi. 

The redshift dependence of Ab/bL, 
(Afo/6^)a[2/jv£QM(fc)g]^^ is shown in Fig. [10] where 
M > W^^Mq/H and scales k < 0.026 h/Mpc were used. In 
applying the correction 5c/(t\i — > q5c/cr\,i to A6/(&g — 1) 
we have actually corrected b{f , i.e. we have employed the 
same approximation used for the green dot-dashed line of 
Fig. 8, giving Eq. Eq.(9) in fact is only the consequence 
of our correction to the Gaussian halo bias. Note that the 
approximation Zf ~ 20 we employed here is expected to 
hold for rare-massive-halos and Fig. 8 shows that this is 
a good approximation. A detailed study of the dependence 



of the non-Gaussian halo bias correction on the formation 
redshift of the halos will be presented elsewhere. 

There seems to be an indication that the g-correction 
factor for the large-scale bias correction may slightly depend 
on the value of /jvz,: in particular the figure shows that it 
could be slightly smaller than q for Jnl large and negative 
and smaller for fNL large and positive. This is not unex- 
pected: the presence of non-Gaussianity may alter the dy- 
namics of non-spherical colla pse (e.g., through tidal forces 
- see e.g., iDesiacaues I (|2008l ) - or by significantly changing 
the redshift for collapse with respect to the Gaussian case). 
At this stage, however, this trend is not highly significant 
and further study will be left to future work. 

Finally, we show the scale dependence of Eq. (|9]), 
(A6/6f )s[2/ivL<5c(2)g]"\ in Fig. [11] The thin lines corre- 
spond to different redshifts and the thick black line to 
their average. The dotted line is the theory prediction with 
q = 0.75. Note that there is an excellent agreement on the 
scales of interest, e.g., k < 0.03 /i/Mpc. On smaller scales 
the effect of non-Gaussianity is very small and the measure- 
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Figure 9. Weak mass dependence of Ab/b at scales k < 0.03 
h/Mpc. Different lines correspond to different redshift snapshots 
between 2 = and z = 1.5. The overall normalization is arbitrary. 
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Figure 10. The redshift dependence of the non-Gaussian correc- 
tion to the halo bias: points are the values measured from the 
simulations, lines are the theoretical predictions, Eq. ||9]l. Only 
k < 0.026 h/Mpc were used. 
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Figure 11. Scale dependence of Eq. (|9}. The thin lines corre- 
spond to different redshifts for halos with mass above IO^^Mq/H 
and the thick black line is their average. The dotted line is the 
theory prediction with q = 0.75. At A; > 0.03h/Mpc the effect 
of non-Gaussianity is very small and the measurement become 
extremely noisy. 



ment become extremely noisy. The se results are in qualita- 
tive agreement with the findings of IPillepich et all (|200^ . 

We conclude that Eq. ((9)l, with q ~ 0.75, provides a 
good fit to non-Gaussian simulations. 



7 COMPARISON WITH PREVIOUS WORK 

After the discussion of §2, it should be clear that if fi^if^ 
were used in the theoretical predictions or S3 were not lin- 
early extrapolated to 2 = 0, then any constraints on non- 
Gaussianity so obtained would have to be re-scaled by a 
factor ~ 1.3. This seems to be the case of some work in the 
literature. On the other hand the q-correction factor effec- 
tively introduces a re-scaling of a factor ~ 0.75^^^^ = 0.65 
for the mass function case and 0.75 for the bias case. It 
is a coincidence that, for the halo bias, 1.3 x 0.75 ~ 1, 
thus the /jvL normalization mistake cancels out with the 
spherical collapse approximation error. This fortuitous can- 
cellation does not happen to the same level in the mass 
function 0.75^''^ x 1.3 ~ 0.8, explaining perhaps some of 
the claimed discrepancy of the simulations with the ana- 
lytic mass function predictions, and the claimed agreement 
with the halo bias predictions. Another possible source of in- 
accuracy would be an inconsistent treatment of the redshft 
evolution of 5c and S3 (see discussion is §2). 

In Fig. [T5] we compare our theoretical predi ctions 
with the results prese nted in IPillepich et al.l l|2008l ) and 
iDesiacgues et af] (I2OO8I ). The left panel shows our simula- 
tion results at 2 = for /nl = ±100 and our theoreti- 
cal predictions . Addi tionally we show the fit presented by 
IPiUepich et~al] (120081 ). Eqs. (8) and (9), evaluated for the 
suitable values of /nl accounting for the different notations 
for /nl- We also adopt our cosmological parameters when 
converting ctm to M. The right panel sho ws the simulation 
results presented in IPillepich et af] (|2008l ) at z=0 and their 
fitting formula at 2 = and 2 = 1. We over plot our the- 
oretical models evaluated for their cosmological parameters 
and for the correspo nding values of /nl- Mo reover we add 
the data points from IDesiacgues et al.l (|2008l ) for 2 = and 
2 = 1, suitably rescaled by the differences of the /nl value 
used. The mass function fits of Fig. 12 differ for large masses, 
in the regime where simulations errors become large; the 
fits are however consistent given the individual points error- 
bars. 

Our theoretical formulae for the non-Gaussian mass 
function (Eqs. 5, 6, 7) and for the non-Gaussian halo bias 
(Eq. 9) are physically motivated expressions that have been 
tested on N-body simulations. They have the advantage 
over fitting formulae that they can be more robustly in- 
terpolated and extrapolated to cosmologies and parame- 
ters that have not been directly simulated and they are 
more robust over parameters ranges where the simulations 
have low signal-to-noise. Compared to simple fitting for- 
mulae, Eq. 6,7 and 9 have the disadvantage that they re- 
quire the calculation of some numerical integrals. To over- 
come this, we supply tabulated values for 5*3^^, am and 
um for a WMAP5 cosmology in the range of interest at 
www.ice.csic.es /personal / verde / nongaussian.html. 

The g-correction we find here has implications for 
previously reported and forecasted constraints on non- 
Gaussianity. In Table 1 we report present and forecasted 
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Figure 12. Left panel: comparison between simulations points and fits from this work and the polynomial fit of iPillepich et aL ] 1I2OO8I) . 
Right panel: Our correction to the Matarrese, Verde & Jimenez (2001) and LoVerde et al (2008) non-Gaussian mass function fits, the 
polynomial fit of Pillepich, Porciani & Zahn (2009) and points from Desjaques et al (2009) (black squares) and Pillepich, Porciani & 
Zahn (2009) (green triangles). See text for more details. 



constraints on /nl from the literature rescaled to /^l 
and corrected for our factor q. 

This confirms that constraints on /jvl achievable 
using the non-Gaussian halo bias are competitive with 
CMB constraints (/jvl ~ 5 for P lanck an d Jn l ~ 3 
for a CMBPol-type mission, 



lYadav. Komatsu fc Wandelt 



Babich fc Z aldarriaga I l|2004l) : 
" (|2007l)). 



8 CONCLUSIONS 

We have considered (1.2 Gpc/h)'^ size and 960^ particles N- 
body simulations with non-Gaussian initial conditions, with 
non-Gaussianity parameter /jvl = ±100, /jvl = ±200 and 
a reference Gaussian simulation (/jvl = 0). The clustering 
properties and the abundance of the simulation's halos were 
then compared with independent simulations and theoreti- 
cal predictions. We find good agreement between different 
simulations, indicating that the initial conditions set-up is 
under control. We find that the Press-Schechter-based de- 
scription of the non-Gaussian correcti on to the Gaussian 
mass function of lMat arrese et all (I2OOGD and iLo Verde et al.l 
l|2008l ) is a good fit to the simulations, provided that: 
a) The Press-Schechter-based description is used to com- 
pute the ratio between Gaussian and non-Gaussian mass 
function. 

h) The critical density 5c is corrected to account for non- 
spherical collapse dynamics. 

This is summarized in our Eq. ((5} an d in Eqs. (|6l) and Q for 
the non-Gaussian mass fu nctions of lMatarrese et al] (|200G| ) 
and ILo Verde et al] (|2008l ). respectively. For large thresholds 
this correction is equivalent to a re-scaling of the spherical 
collapse threshold 5c x ^ where q — 0.75. The q-correction 
is thus equivalent to a reduction of /jvl by a factor ~ 1.5 be- 
cause in the mass function, to leading order /jvl multiplies 

We find that t he non-Gaussian halo bias p r escription of 
iDalal et all l|200'iD . iMatarrese fc Ver"d3 l|2008l ). ISlosar et all 



Table 1. Current and forecasted constraints on f^^^^ 



Data/method 



Measurements 
/nl , 1 — 2(7 errors 



reference 



Photo LRG-bias 
Spectro LRG-bias 
QSO - bias 
combined 



,+54+85 
*-101-331 



93 



11 



+74+139 



-26+47 



37 



+ 23+42 



Slosar et al. 2008 
Slosar et al. 2008 
Slosar et al. 2008 
Slosar et al. 2008 



NVSS-ISW 
NVSS-ISW 



140 



+647+755 



Slosar et al. 2008 



-337-1157 

272 ± 127 (2-cr) Afshordi&ToUcy 2008 



Data/method 



Forecasts 
A/nl(1-'t) 



reference 



BOSS-bias 
ADEPT/Euclid-bias 
PANNStarrs-bias 
LSST-bias 



18 
1.5 
3.5 
0.7 



Carbone et al 2008 
Carbone et al. 2008 
Carbone et al. 2008 
Carbone et al. 2008 



LSST-ISW 



10. 



Afshordi&ToUcy 2008 



(|2008h and lAfshordi fc TollevI | 20081 ) provides a good de- 
scription of the scaling of the large-scale halo clustering of 
the simulations. In particular, we have tested separately the 
predicted redshift, scale and /jvl dependence. The overall 
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amplitude of the effect, however, should be corrected by a 
factor ~ q which can also be understood in the context of 
ellipsoidal collapse or as a modification to the excursion set 
ansatz and the sharp-k space filtering (see Eq. (H}). There 
is an indication that this correction may be slightly depen- 
dent on /nl- This is not unexpected, but the signal-to-noise 
of the effect is too small in the current simulations to draw 
robust conclusions. We also find that on large {k < 0.03 
ft/Mpc) scales, as expected, the fractional correction to the 
non-Gaussian halo bias is independent of mass. On smaller 
scales a dependence on mass is expected, but the simula- 
tions do not have sufficient signal-to-noise to verify it. The 
g-correction to the non-Gaussian halo bias modifies current 
and forecasted constraints reported in the literature as in- 
dicated in our Table 1. 

The formulae we presented here for the non- 
Gaussian mass function (Eq. 5, 6,7) and non-Gaussian 
halo bias (Eq. 9) are physically motivated expressions 
which provide good fits to a suite of N-body simu- 
lations. As such, they can be more robustly interpo- 
lated and extrapolated than simple fitting functions (in 
www.ice.csic.es/personal/verde/nongaussian.html we pro- 
vide useful quantities for ease of use of these equations). 
We confirm that the non-Gaussian halo bias offers a robust 
and highly competitive test of primordial non-Gaussianity. 
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